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ABSTRACT 

During accretion a neutron star (NS) is spun up as angular momentum is transported through its surface 
layers. We study the resulting differentially rotating profile, focusing on the impact this has for type I X-ray 
bursts. The predominant viscosity is likely provided by the Tayler-Spruit dynamo, where toroidal magnetic 
field growth and Tayler instabilities balance to support a steady-state magnetic field. The radial and azimuthal 
components have strengths of ~ 10 s G and ~ 10 10 G, respectively. This field provides a Maxwell stress on the 
shearing surface layers, which leads to nearly uniform rotation at the depths of interest for X-ray bursts (near 
densities of w 10 6 g cm" 3 ). A remaining small shear transmits the accreted angular momentum inward to the 
NS interior. Though this shear gives little viscous heating, it can trigger turbulent mixing. Detailed simulations 
will be required to fully understand the consequences of mixing, but our models illustrate some general features. 
Mixing has the greatest impact when the buoyancy at the compositional discontinuity between accreted matter 
and ashes is overcome. This occurs at high accretion rates, at low spin frequencies (when the spin is small, the 
relative speed of the accreted material is larger), or may depend on the ashes from the previous burst. We then 
find two new regimes of burning. The first is ignition in a layer containing a mixture of heavier elements from 
the ashes. If ignition occurs at the base of the mixed layer, recurrence times as short as ~ 5 - 30 minutes are 
possible. This may explain the short recurrence time of some bursts, but incomplete burning is still needed to 
explain these bursts' energetics. When mixing is sufficiently strong, a second regime is found where accreted 
helium mixes deep enough to burn stably, quenching X-ray bursts. We speculate that the observed change in 
X-ray burst properties near one-tenth the Eddington accretion rate is from this mechanism. The carbon-rich 
material produced by stable helium burning would be important for triggering and fueling superbursts. 
Subject headings: accretion, accretion disks — stars: magnetic fields — stars: neutron — X-rays: bursts — 
X-rays: stars 



1. INTRODUCTION 

As neutron stars (NSs) in low mass X-ray binaries accrete 
material from their companions they are expected to be spun 
up by this addition of angular momentum, possibly becomin g 
millisecond pulsars (Bh attacharya & van den Heuvell [T991). 
This suspicion has received suppor t by the discovery of ac- 
cretio n driven millisecond pulsars ( Wiinan ds & van der Klisl 
1998), as well a s the millisecond oscilla tions seen during type 
I X-ray bursts dChakrabartv et alj|2003l) . t he unstable ignition 
of the accumulating fuel (for reviews, see iLewin et alJll995b 
IStrohmayer & Bildstenll2006l:lGailoway et alJl2006h rThe ma- 
jority of the rotational kinetic energy of the accreted material 



is dissipated at low densities (~ 10 
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boundary layer (for example, llnogamov & Sunvae\ill999h . 
Nevertheless, angular momentum must be transported into the 
NS interior if it is to be spun up, even if at times th e angular 
mome ntum could be radiated as gravitational waves (Bildsten 
1998b). This transport implies a non-zero, albeit small, shear 
throughout the outer liquid parts of the NS. Such shearing 
may lead to viscous heating as well as chemical mixing at 
depths far below the low density boundary layer where the 
majority of the shearing occurs. 

Such a pi cture was previously investigated by Fuiimoto 
d!988l 119931) for both accreting NSs and white dwarfs. His 
studies demonstrated that the shearing profile of a NS could 
be large enough that mixing through baroclinic instabilities 
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may be important at the depths of type I X-ray bursts. Such 
a result is compelling because it may help to explain some of 
the remaining discrepancies between the theory and observa- 
tions of these bursts. Though a simple limit cycle picture has 
been successful in qualitatively explaining the bursts' primary 
characteristics, including their energies (10 39 - 10 40 ergs), re- 
currence times (hours to days), and d urations (~ 10- 100 s) 
dFuiimoto et al Jl 1 98 It iBildstenll 1 998al). outstanding problems 
still remain (IFujimoto et all Il987t Ivan Paradijset all 11988: 
lBildstedl2000h . Chief among these is the critical accretion 
rate required for stable accumulation (which we discuss in 
more detail in $6]). Another problem is the occurrence of suc- 
cessive X-ray bursts with recur rence times as short as ~ 10 
minutes (Gallo way et all [20061 and references therein), for 
which there has been sp eculation that this c ould be due to 
mixing during accretion (Fuii moto et alJfl987l) . 

In this present work we reassess the importance of angular 
momentum transport and the resulting mixi ng. We find tha t 
the hydrodynamic instabilities studied by Fuii moto! d 19931) 
are insufficient to prevent strong shearing of magnetic fields. 
This leads to generation of the Tayler-Spruit dynamo d Spruit! 
|2002|) . where toroidal field growth is balanced by Tayler in- 
stabilities to create a steady-state magnetic field, which pro- 
vides a torque on the shearing layers that is larger than any 
purely hydrodynamic mechanism. The result is nearly uni- 
form rotation and little viscous heating, too little to affect 
either X-ray bursts or superburst s (thermonuclear ignition of 
ashes from previous X-ray bu rsts,[Cumming & Bildste rj2001l : 
IStrohmayer & Brownll2002l) . Turbulent mixing is found to 
be non-negligible and in some cases may mix fresh mate- 
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rial with the ashes of previous bursts. The key is whether 
the strong buoyancy due to the larger density of the ashes be- 
low can be overcome. When mixing occurs, we find two new 
regimes of burning for an accreting NS. The first is that fuel 
can be mixed down to depths necessary for premature unsta- 
ble ignition. The timescale for ignition of such bursts is short 
enoug h (~ 5-30 minutes) to explain the short recurrence time 
bursts dGallowav et al.l 120061) . The second is that if material 
is mixed to sufficiently large depths (and therefore tempera- 
tures) it can burn stably, ceasing X-ray bursts altogether. Such 
a mechanism may be responsible for the quenching of X-ray 
bursts seen at surprisingly low ac cretion rates for many atoll 
class NSs (ICornelisse et al.ll2003l) . We derive an analytic for- 
mula that estimates what spins and accretion rates necessary 
to transfer between these two regimes (eq. IfTTl ). 

1.1. The Basics of Angular Momentum Transport 

Before beginning more detailed analysis, it is useful to 
present some general equations for angular momentum trans- 
port. This demonstrates why we expect the NS interior to be 
shearing and explicitly relates the torque from the accreting 
matter to this shear. Figure [T] schematically shows the ex- 
pected rotation profile from the accretion disk, through the 
boundary layer, and into the NS interior. Material accreted 
at a rate M reaches the NS surface with a nearly Keplerian 
spin frequency of Q K = (GM/R 3 )^ 2 = 1.4 x 10 4 s" 1 m\^R~ V2 , 
where M 14 =M/\AM Q and R b =R/\Q b cm, which has a ki- 
netic energy per nucleon of w 200 MeV nuc" 1 . The major- 
ity of this energy is dissipate d in a boundary layer of th ick- 
ness H B l <C R (as studied by Inoga mov & Sunyaevl 19991) and 
never reaches far into the NS surface. Nevertheless, angular 
momentum is being being added at a rate of M(GMR) 1 / 2 , so 
that a torque of this magnitude must be communicated into 
the NS. This implies a non-zero shear rate in the interior liq- 
uid layers, down to the solid crust. In the present work we are 
interested in the shear at the depths where X-ray bursts ignite, 
near p w 10 6 g cm" 3 , which is well below the boundary layer. 
Though we focus on a magnetic mechanism for angular mo- 
mentum transport, we expect the boundary layer and densities 
up to « 5 x 10 3 g cm" 3 to be dominated by hydrodynamic in- 
stabilities (which we discuss in more detail in ® . 

The pressure scale height at the depth of X-ray bursts is 
H m 30 cm <C R, which allows us to use a Cartesian coordi- 
nate system with z as the radial coordinate. Thi s is far below 
the depths considered by Inogamov & Sunyaev (1999) when 
they investigated the uneven covering of the NS surface by 
accreted fuel. Furthermore, all of the transport mechanisms 
we consider work most efficiently in directions perpendic- 
ular to gravity (because no work is performed), thus it is a 
good approximation to consider the surface layers as concen- 
tric spheres, each with constant Q. Transfer of angular mo- 
mentum is redu ced to a one-dimensional diffusion equation 
dFujimotoll 19931) 
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where fl is the NS spin frequency and v is the viscosity. The 
total time derivative is given by d/dt = d/dt + V^d/dz, where 
Vadv is the advecting velocity of the fluid in an Eulerian frame. 
In steady-state, we take d/dt = and set V a d v = -M/{AnR 2 p) 
so that 
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FIG. 1. — A diagram demonstrating the rotation profile of material from 
the accretion disk, through the boundary layer, and into the NS. Material 
reaches the NS surface with a Keplerian spin frequency, Ck, the majority of 
which is dissipated within the boundary layer (gray region) with thickness 
//bl <C R- Nevertheless, angular momentum is still being added to the NS, 
and a torque must be communicated through the NS. This implies a small, 
non-zero, amount of shear throughout the NS interior. 

Integrating from the surface where the spin is J7k (due to the 
disk) down to a depth z where the local spin is f2, and assum- 
ing dVl jdz = at the surface, 

-MR 2 tt K +MR 2 n = -4npvR 4 dn/dz. (3) 
Taking the limit -C we find, 



4irR 3 pvqn = MR 2 n K , 



(4) 



M d ~ 1 d ( Adti 
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where q = d\og£l/d\ogz is the shear rate. This equation 
shows that q > when angular momentum is transported in- 
ward. In general, we will find that q is rather small (< 1) at 
the depths of interest. Nevertheless q is large enough to acti- 
vate instabilities that help to transport angular momentum, as 
well as mix material. 

1.2. Outline of Paper 

We begin by comparing and contrasting some well-known 
hydrodynamic instabilities in ^2] and estimate the resulting 
shear rates. In |3j we discuss the consequences that this shear- 
ing has on a magnetic field, which motivates implementation 
of the Tayler-Spruit dynamo. In we calculate accumulat- 
ing NS envelopes without the effects of viscous angular mo- 
mentum transport. This enables us to judge when such effects 
must be incorporated. We calculate models including mixing 
in $5] We conclude in $6]with a summary of our results and a 
discussion of type I X-ray burst and superburst observations. 

2. HYDRODYNAMIC VISCOSITY MECHANISMS 

In the following sections we discuss hydrodynamic insta- 
bili ties. This is not a n exhaustive survey (for further details, 
see lHeger et al.l l2000). but rather is meant to summarize those 
instabilities that are most crucial to our problem, so as to set 
the context for the magnetic transport mechanism we study 
later. 

2.1. Kelvin-Helmholtz Instability 

The Kelvin-Helmholtz instability (also referred to as the dy- 
namical shear instability) is governed by the Richardson num- 
ber 

N 2 
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where N is the Brunt-Vaisala frequency, which is composed 
of contributions from thermal and compositional buoyancy, 



n 2 =n$+nI. 



The thermal buoyancy is given by 



N 2 = l^ 



V ad - 



fd\nT\ 



(6) 



(7) 



where g = GM/R 2 = 1.87 x 10 14 cm %- 2 M\ A Rf is the surface 
gravitational acceleration (ignoring the effects of general rela- 
tivity), xq = dlnP/dlnQ, with all other intensive variables set 
constant, V a( j = (<91nr/<91nP) ac j is the adiabatic temperature 
gradient, the asterisk refers to derivatives of the envelope's 
profile, and the pressure scale height is, 



H = P/pg = 33.1 cm/iijjTg, 



(8) 



where 7s = T/IO 8 K. For these analytic estimates we assume 
an ideal gas equation of state and use a pure helium composi- 
tion with a mean molecular weight of p = 1.33//1.33. We omit 
the scalings with mass and rad ius to simplify presentation . 
The compositional buoyancy is dBildsten & Cummingll 19981) 



N 2 = - 



d\np e 
dlnP 



( dlnp/\ 



(9) 



where p e (/x,) is the mean molecular weight per electron (ion). 
Linear analysis shows that Kelvin-Helmholtz instability oc- 
curs when Ri < 1/4, which develops into strong turbulence 
that readily transports angular momentum. This result as- 
sumes that thermal diffusion can be ignored for the unstable 
fluid perturbations, in other words, that the perturbations are 
adiabatic. 

Fluid perturbations with a characteristic size L and speed V 
become non-adiabatic when the timescale for thermal diffu- 
sion, L 2 /K, where K is the thermal diffusivity, is less than 
the timescale of the perturbation, L/V. The rati o of these 
two t imescales is the Peclet number, Pe = VL/K dTownsendl 
1958). The restoring force provided by thermal buoyancy 
is weakened when Pe < 1, which requires the substitution 
of Nj — > PeNj and promotes instability for regions where 
Nt > N^. Thermal diffusion is most efficient at small length- 
scales, which motivates setting LV/vk = Re c dZahnl ll992). 
where Vk is the kinematic viscosity and Re c is the critical 
Reynolds number for turbulence, which is of order 1000. This 
gives the Peclet number approximately related to the Prandtl 
number, Pr, by Pe ks Re c Pr. The tu rbulent perturbations are 
thus non-adiabatic when (Zahn 19921) 



K > v k Re c 



(10) 



In the non-degenerate surface layers the kinematic vi scosity 
is dominated by ions, and has a value of (Spi tzerll 1962b 



1.4 x 10" 3 cm 2 s^p^T^ 2 , 



(11) 



where p(, = p/10 6 g cm 3 , and we assume a Coulomb loga- 
rithm of In A = 20. Setting K= 16os,bT 3 /(3c p Kp 2 ), where ctsb 
the Stefan-Boltzmann constant, c p the specific heat, and k the 
opacity, the thermal diffusivity is 

A" = 48.8 cm 2 s"Vi Wm^ 3 , (12) 
where we approximate c p = 5k^/2pm p and scale the opac- 
ity to Ko.04 = ft/0.04 cm 2 g" 1 (the opacity is largely given by 
electron scattering, but is decreas ed by degeneracy effects, see 
lPaczvnskifl 983; Bil dsten|[T998ah . Substituting equations O 



and (fT2l into equation ( [Tol l, we find that the perturbations are 

non-adiabatic at depths of p < 3 x 10 7 g cm" 3 T^ 1 . The new 
"secular" Richardson number associated with this limit is, 



Rh 



K q 2 il 2 ' 



(13) 



When Ri s < 1/4, the so-called "secular shear instability" 
arises. 

The competing effects of accretion increasing q versus tur- 
bulence developing when Ri s < 1/4 (and decreasing q) drive 
the surface layers toward marginally satisfying Ri s =1/4 (as- 
suming at this moment that the sole viscous mechanism is the 
Kelvin-Helmholtz instability). This expectation is bo rne out 
in the white dwarf studies of I Yoon & Langer (2004). Thus 
we can trivially esti mate the q due t o this mechanism. The 
thermal buoyancy is (Bildsten 1998a) 



1/2 



= 9.2x 10 5 s 



-1,1/2 -1/2 



(14) 



We substitute Ri s =1/4 into equation dT3b . and assuming 
Re c = 1000, solve for a shear rate of 



(15) 



OOQ 1/2 l/2„-3/4 - 
4KH = 223 KqWV r 8 ty).li 

where fio.i = f2/0.1f2K- A shear rate this large would promote 
prodigious viscous heating, as well as ample mixing, but as 
we soon show, such a large shear is prohibited by other insta- 
bilities. 

2.2. Baroclinic Instability 

Another important hydrodynamic instability that has been 
studied extensivel y for accreting d egenerate stars is the baro- 
clinic instability (Fuiimotol ll993h . which we qui ckly sum- 
marize here. The interested reader should consult Fuiimoto 
(119871119881) for further d etails (also see the discussion in 
ICumm ing & B ildstenll2000l) . 

The baroclinic instability arises because surfaces of con- 
stant pressure and density no longer coincide if hydrostatic 
balance is to be maintained when differential rotation is 
present. In such a configuration, fluid perturbations along 
nearly horizontal directions are unstable, though with a suffi- 
cient radial component to allow mixing of angular momentum 
and material. The instability can roughly be broken into two 
limits, dependin g on a critical baroclinic Richardson number 
dFujimotolll987h . 



A'/ B( =4 1 ^] 



H 



N 



: 8.5xl0 3 M i.337F 1 ^o.i- (16) 



When Ri > RiBC, Coriolis effects limit the horizontal scale of 
perturbations. This results in tw o parametrization s for viscos- 
ity estimated from linear theory (Fuiimoto 1993), 



1 1 




H 2 n, Ri<Ri BC , 



(17) 



H 2 n, Ri>Ri 



BC 



3 RiV 2 

where a factor of order unity is usually included in these pre- 
scriptions, called chbc, to account for uncertainty in how lin- 
ear theory relates to the saturated amplitudes of the instability. 
For simplicity, we set asc = 1 in our analysis below. 

By substituting v^c into the angular momentum equation 
(eq. ID), we solve for the shearing profile. Since we are 
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interested in how the shear relates to the bursting properties, 
it useful to write these results in terms of the local accretion 
rate m = M/{4irR 2 ), which is typically parametrized in terms 
of the local Eddington rate 



'«Edd = 



2m p c 



(l+X)Ra Th 
1.5 x 10 s g cm" 2 s" 1 
(l+X)R 6 



(18) 



where X is the hydrogen mass fraction and aTh is the Thomson 
cross-section. The Richardson number in each case is then 



Ri 



2.2 x 10 3 ^ipeT^m^ Q .i, Ri < Mbc 
. 4.4 x 10 3 ^pf^^Ml » > Kec 



(19) 

where mo. l = ni/O.lmrjdd an d we h ave use d the hydrogen 
deficient (X = 0) value for «?Edd (eq. fT8"1l ). The transi- 
tion to the case Ri > Rise occurs roughly at depths p > 

4 x 10 6 g cm" 3 T 8 5,/2 . The shear rate for each case is 



<7bc = 



(20) 



This demonstrates that generally q B c <C ^kh> so that the baro- 
clinic instability triggers before the Kelvin-Helmholtz insta- 
bility. This prevents the shear rate from ever becoming large 
enough for the Kelvin-Helmholtz instability to operate at 
depths of p > 6 x 10 4 g cm" 3 . 

2.3. Other Hydrodynamic Instabilities 

In addition to the two instabilities described above, 
there are a number of other possibilities incl uding, but 
not li mi ted to, Eddington-Sweet circulation dvon Zeipel 
ll924alrjtlBaker & Kippeniiiin|[T959l). Solberg-H0iland insta - 
bility iwasiutvnskil 1 946t iTassoulll 1 978t| Endal & Sofia 19781) 
Goldreich- S chuber t-Fricke instability (Goldreich & Schubert 
1 1 9671: iFrickell 19681) and Ekman pumping (Pedlos ky|1987l) . At 
this time we avoid assessing each of these individually. As we 
show below, magnetic fields are likely to dominate and are in- 
teresting since they have received the least attention in past 
works. 

3. THE IMPORTANCE OF MAGNETIC EFFECTS 

Given the larger than order unity shear rates derived above, 
we estimate the consequences this has for a magnetic field. 
The point we want to emphasize is that even a reasonably 
small field will be wrapped by the shear flow until it becomes 
dynamically important. 

Assuming shellular rotation, a component of radial field is 
stretched to have a toroidal component, = nB r , where n is 
the number of differential turns, given as n = qflt, and t is 
the duration of the shearing. The toroidal field growth is very 
fast. For the Kelvin-Helmholtz case, ~ B r in merely ~ 
10~ 6 s (~ 10~ 4 s for the baroclinic case). As the toroidal field 
becomes larger it exerts an azimuthal stress on the shearing 
layer equal to B, B^ /Air, which can be written as an effective 
viscosity, v e , 



B r B^ 

47T 



= pv e qfl. 



(21) 



In a timescale t m H^/v e this torquing significantly decreases 
the spin of the layer, where = dz/dlnfl = R/qis the shear- 
ing scale height. Setting B^ = qiltB n we solve for the critical 
initial radial field needed to affect the shearing, 



Br.crit = (47rp) 



1/2 



He, 



: (4irp) 



1/2 



R 

7r 



(22) 



so that t is basically the Alfven travel time through a shearing 
scale height. Using p 10 6 g cm" 3 as burning density and t w 
1 hr, a fiducial timescale for accumulation, our estimates for 
<7kh and q^c imply B naA ~ 10 4 G and ~ 10 s G, respectively 
Any initial field larger than this gets so wound up that Lorentz 
forces alter the shearing profile. 

It is possible that the intrinsic magnetic field of the NS may 
be large enough that mag netic stresses n ever allow the shear 
rates to become so large. ISpruij dl999l) argued that this de- 
pends on the "rotational smoothing time," the timescale for 
non-axisymmetric components of the magnetic field to be ex- 
pelled as differential rotational brings field with opposite po- 
larities together. This timescale is estimated as 



f3R 2 ir 2 



\r,Vl 2 q 2 



1/3 



(23) 



where 77 is the magnetic diffusivity, w hich in the non- 
degenerate limit is given by dSpitzerll 19621) 



77 w 0.7 cm s To . 



(24) 



If to is larger than the Alfven travel time through the layer 
(given by t in eq. 1221 ). then there is sufficient time for mag- 
netic torques to act and the shearing is merely a perturbation 
on the magnetic field that is quickly damped away. On the 
other hand, if fn is sufficiently small then the shear dominates 
and the magnetic field is made axisymmetric on a timescale 
fa. Substituting fn into equation ( 1221 . we estimate the critical 
field, below which shear dominates, 



B Q = (4irp) 



1/2 



TJ&R 



1/3 



-1/3 



3tt 2 



(25) 



The lack of persistent pulsations fro m accreting NSs implie s 
a dipole field strength < 5 x 10 7 G dPiro & Bildstenll2005h . 
We therefore consider it plausible that they may have intrin- 
sic magnetic fields < Bq. In this case, we expect a very 
axisymmetric magnetic field to be created in a timescale of 
tfi ~ 200 s, which is then wrapped until it becomes dynami- 
cally important. 

3.1. The Tayler-Spruit Dynamo 

As the toroidal magnetic field continues to wrap, it becomes 
increasingly important to the dynamics of the shearing and is 
also subject to magnetohydrodynamic instabilities. The com- 
bination of these effect s have been shown to give rise to the 
Tayler-Spruit dynamo (Spruit 2002). In this picture, shear- 
ing grows the toroidal field, which then initiates Tayler in- 
stabilities (non-a xisymmetric, pinch-like in stabilities includ- 
ing stratific ation. iTaylerll 1 973"t ISpruitll 1 9991) . This turbulently 
creates poloidal field components that once again shear to be 
toroidal. This cycle continues, creating a steady-state field. 

The minimum shear needed for this process to operate can 
be argued simply. (See Spruit 2006 for a mathematical deriva- 
tion that uses the dispersion relation from Acheson 1978.) We 
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note that iDenissenkov & Pinsonnea ult ( 2007]) have recently 
given an alternate prescription for this same mechanism using 
solely heuristic arguments. Since their results have not been 
shown consistent with a more rigorous mathematical analysis 
we consider Spruit's conclusions more reliable at this time. A 
vertical perturbation, l z , is limited by buoyancy forces to be 
(eq. [6] from Spruit 2002) 



U < Ruja/N, 



(26) 



where cja = B/[(Ai:p) x l 2 K\ is the Alfven frequency. At small 
lengthscales, magnetic diffusion damps out perturbations. In 
the limit of ^> wa, which we are considering, the Tayler 
instability growth rate is as = a; A /fi, so that 



I 2 > r)/(TB = rjVt/uj 2 h . 



(27) 



Combining these two relations gives the minimum u>a needed 
for the dynamo to act, 



o; A \ 4 V (N 



(28) 



During the timescale for Tayler instability, a B l ,B r is stretched 
into Z?0 by an amount 



B 4> = cr B l qflB r . 



(29) 



The largest amplification is achieved for magnetic fields that 
extend the largest radial lengthscale available, so that assum- 
ing equation (T2"6l ) is marginally satisfied along with the induc- 
tion equation we find B r jB^ = U/R = loa/N. Combining this 
with equation (|29l we obtain q = (N /ST)(uja/ Q)- Substituting 
this into equation (1281 1, we find (Spruit 2002) 



n ) \r 2 n) 



1/4 



(30) 



as the minimum shear needed for the Tayler instability to op- 
erate. T hough this result is consistent with more rigorous 
analysis dSpruitl2006l) . it should be viewed with some caution 
as we apply it to accreting NSs. The thin shell geometry we 
consider is quite different than the spherical geometries typi- 
cally us ed when invoking the Tayler-Spruit dynamo. Simula- 
tions by Braithwaite (2006) demonstrate that the Tayler insta- 
bility is strongest along the rotation axis, which is only real- 
ized at the poles in the NS case. Since we find the dynamo to 
be so much stronger than any hydrodynamic transport mech- 
anisms, we consider it to be a reasonable approximation for 
these magnetic effects, even if its strength is decreased due to 
geometry. 

The value of N used to evaluate equation d30l depend s on 
what is supplying the buoyancy. We follow Spruit] d2002l) and 
separately consider cases of 3> Nt and <CWj-, denoted 
as case and 1, respectively. In case 1 non-adiabatic effects 
become important when rj/K < 1, and we must take N T — > 
{rj/K)Nj (analogous to the above analysis of the secular shear 
instability). For this case, which dominates for most of the 
envelope, we find 



-o io k /4 3/ V 9/2 rr 7/4 



(31) 



Both <7kh and q^c are considerably above this, thus the Tayler- 
Spruit dynamo activates long before the onset of purely hy- 
drodynamic instabilities. 

By a ssum ing that equation (|28| | is marginally satisfied 
Spruit (2002) derived the steady-state field strengths. We 



summarize the relevant prescriptions needed for our study. 
The steady-state azimuthal and radial field components are 

B^^irp^Rqn 2 /^, (32) 

2 



B,o = q — fi^o, 



for case 0. For case 1, 



fi i=(4^o) 1/2 RfV /2 



L/8 



K 



n 



Nt) \R 2 N t 



I- 



1/8 



B r \ = 



0i. 



S N T J \R 2 N T 
The effective viscosities, as defined by equation ( 1211 ). are 

n 



u =R 2 nq 2 



n\ 1/2 / k 



1/2 



(33) 

(34) 
(35) 

(36) 
(37) 



^N T J \R 2 N t , 

Although these viscosities are appropriate for angular mo- 
mentum transport, mixing of material is less efficient since it 
requires expending work to exchange fluid elements (versus 
just exerting shear stresses). The mixing diffusivities are 

6 



D x =R 2 n q 



-)"'( 

Nt) \R 2 N t 



K 



3/4 



(38) 



(39) 



which is just equal to the effective turbulent magnetic diffu- 
sivity. In the Appendix we show that these prescriptions are 
consistent with energy conservation considerations. 

3.2. Shearing Profile 

The Tayler-Spruit dynamo transports angular momentum, 
causes viscous heating, and mixes material. We follow the 
procedure we used for the baroclinic case and assume steady- 
state angular momentum transport. This is a good approxima- 
tion when the viscous timescale, f v j sc = H 2 /u, is less than the 
timescale of accretion. For these initial estimates we focus 
on the Nt S> iV M limit (case 1 ) since this dominates except at 
compositional boundaries (which we revisit in ©. Substitut- 
ing V\ into equation (01, the shear rate for the Tayler-Spruit 
dynamo is 



m 17k 

~Rp~n~n 



1/2 



R 2 N 



1/2 



Scaling to values appropriate for accreting NSs, 



<yrs = 0.38 K '^Tg 



2 • o" 5 / 2 

'«o.i"ni • 



(40) 



(41) 



Note the scalings with m and fi. At high m angular momen- 
tum is fed into the star faster, creating more shear. At low 
fl the shear is greater because of the larger relative angular 
speed between the accreted material and the NS. These are 
generic feature we expect for any viscous mechanism (com- 
pare the scaling of gxs with ^kh and ^bc from eqs. fTBI 
and l20l . respectively). The viscosity that gives the small- 
est shear rate will likely be the most important at a given 
depth. Using this criterion, we find that the Tayler-Spruit dy- 
namo is dominant for densities p > 3 g cm" 3 T g 5//2 »ig { , 
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or p > 5 x 10 3 g cm" 3 for T w 5 x 10 6 K. At shallower 
depths, Kelvin-Helmholtz instabilities damp the shear, which 
is consistent with the use of hydrodynamic instabilities by 
Inoga mov & Sunvaevl ([1999) for understanding the initial 
spreading of accreted material in the boundary layer. 
The steady-state magnetic field components are 

B, = l.3 x 10 10 G /^VW^Xf: (42) 

B r = 2.l x 10 5 G K^pfX^oXu- (43) 

ICumming & Bildstenl ( 120001) argue that if B r > 10 6 G it would 
become dynamically i mportant to determ ining the drift of X- 
ray burst oscillations (|Muno et al.ll2002l) . The interaction of 
such fields with the shearing of an X-ray burst from an accret- 
ing m illisecond pulsar has been explored by lLovelace et al.l 
(2007). At high accretion rates, B r increases and comes close 
to this limit, suggesting that these magnetic fields may be 
important for understanding the dynamics of X-ray burst os- 
cillations, even from NSs that do not show persistent pulsa- 
tions (which are preferentially seen at high accretion rates, 
iMuno et alJl2004h . 

3.3. Viscous Heating 

Viscous shearing heats the surface layer, which can also 
be thought of as the rate of magnetic energy destruc- 
tion as the dynamo bu ilds and destroys magnetic field 
(Mae der & Meynel2004l) . The heating rate per unit mass is 

e= l -v{qU)\ (44) 
which for the Tayler-Spruit dynamo becomes 

c ti irvlO -l -l 1/2 -\rr-2 -2 (-1-3/2 

e TS =5.6 x 10 1U ergs g 's K ' M p 6 L T s A Sl { . 

(45) 

To put this number into perspective we express it in terms of 
the energy released per accreted nucleon. This is found by 
multiplying the above result by the total mass per unit area 
down to depth of interest, y m Hp (the column depth), and 
dividing by m. We then find 

=0.13 keV nuc'VTjsS^^'^o.i^o 3 / 2 - ( 46 ) 

In comparison, the thermal energy per nucleon at 10 s K is 
w 10 keV nuc" 1 and burning of helium into carbon releases 
« 0.6 MeV nuc" 1 . We conclude that viscous processes do not 
heat the layer sufficiently to alter X-ray bursts. 

3.4. Turbulent Mixing 

It is not immediately clear how much of an observational 
impact is provided by the magnetic fields and viscous heating 
estimated above. In contrast, as we shall show, shear mixing 
produced by the dynamo could have important ramifications 
for the structure and composition of the surface layers. For 
this reason, we devote most of the remainder of our study to 
investigating the consequences of mixing. 

We first consider some estimates that highlight mixing's im- 
portance. For this, we parametrize the mixing diffusivity as 
«tsA where D is given by equations d38l or ( |39l and «ts is 
a factor of order unity that accounts for uncertainties in the 
Tayler-Spruit prescription. The features we find are general 
enough that other viscosities can be incorporated by increas- 
ing or decreasing a T s. 



The material becomes fully mixed over a scale height, H, in 
a time f m ; x = H 2 / {ctjsD), which gives for case 1 (Nt 3> A^), 

fmix = 4.3 x 10 2 s ajlp^nl^y^Tsm^nl^, (47) 

Accretion is also advecting material downward, which hap- 
pens in a timescale f acc =y//nw pH/rh, 

r acc = 2.2 x 10 3 s MLsaPeTgmo 1 ! . (48) 
The ratio of these timescales is 

tmix „-,„ -l -l 1/4 l/2 3/4 ,.„, 

— = 0.19 a TS p l 3i K ' M p 6 ' Q ' A . (49) 

face 

The scaling with p shows that mixing is generally more im- 
portant than advection at shallow depths (p < 3 x 10 7 g cm -3 ), 
while the scaling with O shows that mixing is more important 
for slower spinning NSs (as expected from our discussion of 
shear rates in 93.21 >. 

4. ACCUMULATING, NON-MIXED MODELS 

We now use numerical calculations to consider the angu- 
lar momentum transport through the surface layers. In this 
section we calculate accumulating models without directly in- 
corporating the effects from mixing. This verifies many of the 
analytic estimates derived above and motivates when mixing 
must be included (which is done for models in 

4. 1 . Shear Profile Calculations 

We calculate the envelope profile for helium accumulating 
on an iron ocean (which represents the iron-peak ashes from 
previous X-ray bursts). As discussed in §12.21 we approxi- 
mate the surface as having a constant gravitational accelera- 
tion, and plane-parallel geometry. In addition to using z as our 
radial coordinate, we find it useful to use the column depth, 
y, defined as dy = -pdz, giving a pressure P = gy from hy- 
drostatic balance. We solve for p using the analytic equa- 
tion of state from Paczyriski (1983). For the liquid phase, 
when 1 < T < 173 where V = (4tt«,/3) 1 ^Z 2 e 2 /(k B T), Z is the 
charge per io n, and w,- is the ion density, w e include the ionic 
free energy of lChabrier & Potekhinl d!998l) . 

During the accumulating phase, a negligible amount of 
helium burning takes place, so we assume that the flux 
is constant and set by heating from the crust. Previ- 
ous studies of accreting NSs have shown that the in- 
terior thermal balance is set by electron captures, neu- 
tron emis sions, and pycnonuclear r eaction in the in- 
ner crust ( Miralda-Escu de et al.l Il990t IZdunik et al.l 119921: 
Bildst en & Brownl 119971: iBrown & Bildstenl 119981: Brownl 
20001 12004 which release » !MeV/m„ » 10 18 ergs g" 1 
(lHaensel & Zduiul fl990l l2003h . Depending on the accre- 
tion rate and thermal structure of the crust, this energy is ei- 
ther conducted into the core or released into the ocean such 
that for an Eddington accretion rate up to « 92% of the en - 
ergy is lost to the core and exits as neutrinos (lBrownll2000l) . 
We therefore set the heating to 150 keV nuc" 1 , giving a flux 
2.2 x 10 21 ergs cm" 2 s" 1 (m)o.i, where (m)o.i is the time aver- 
age accretion rate in units of O.liriEdd and the average is over 
timescales of order the thermal time of the crust (many years). 
For simplic i ty we assume (m) = m. Recent calculations by 
Gupta etail d2006l) suggest that heating is stronger than previ- 
ously thought, but not sufficiently high to qualitatively change 
our results. We ignore the additional flux from compressional 
heating because it only contributes ~ c p T ~ 10 keV nuc" 1 
dUldsten 19981. 
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In a one-zone estimate ignition occurs at the base of the 
helium layer when 



(50) 



dT dT 1 

dFujimoto et al.l 1 198 lb where e^ a is the heating rate from 
triple-a reactions (for which we use the rate from Fushiki & 
Lamb 1987 for numerical calculations), 



^cool ' 



4<T SB T 4 

3ny 2 



(51) 



and the derivatives are both taken at constant pressure. We 
consider models where the column of accreted helium y acc that 
is just at this ignition point. This allows the maximum amount 
of time for angular momentum redistribution during the accu- 
mulating phase. The He/Fe boundary is assumed sharp, since 
the timescale for diffusion between the layers is much longer 
than the timescale of accumulation dBrown et al.ll2002l) . 

We solve for the temperature profile by integrating the ra- 
diative diffusion equation, 



F = 



16<t sb T 3 dT 



(52) 



3k dy 

The opacity is set using electron-scattering ( Paczvnskil 1983b . 
free-free, and conductive opacities (Sch atz et al.ll999l) . There 
is a sharp change in opacity at the He/Fe boundary (because 
the high Z of the iron makes it more opaque), which means 
that the top of the iron layer is convectively unstable. The 
characteristic convective velocity estimated from F « pV^ om 
is Vconv ~ 10 5 cm s" 1 , which is much less than the local sound 
speed of ~ 10 8 cm s" 1 . This means that the convection is very 
efficient, and we therefore simply set {d\nT / d\nP)* = V ac j 
in the convective region. This convective region has not 
been noted in previous studies, though a super-adiabatic tem- 
perature gradient i s ap parent in the fig ures shown in both 
iBrown et all (120021) and lCumm ing (200j). It is not yet clear 
how this impacts the bursting properties of NSs, and we delay 
exploring this in detail for a future study. Nevertheless, we 
must include convection so as to have accurate differentially 
rotating profiles. 

In Figure|2]we plot the temperature profile, K and r), and N 
for three accumulating mod els. The magnetic di ffusivity is set 
using the conductivity from lSchatz et al.1 (I1999I) . The base of 
the accumulating helium is set where unstable helium ignition 
occurs, as designated by the thick long-dashed line in the top 
panel (eq. Il50l ). The convective zone begins just below this 
and is bracketed by solid circles at its top and bottom. The 
convective zone is also seen in the plot of N since Nt is effec- 
tively zero in this region. The change in composition at the 
He/Fe boundary gives a large buoyancy contribution, which 
we estimate as 



f 2. \ '/ 2 



(53) 



where Aln/i w 0.44 is the logarithmic change in the mean 
molecular weight at the boundary. We denote this by a solid 
square for each model in the bottom panel. For the majority of 
the profile the Tayler-Spruit dynamo is given by case 1 (Nt ^> 
N^). Since K ^> 77, the perturbations creating the dynamo are 



non-adiabatic for case 1, as we described in §13.11 The only 
place case (A^ 3> Nt) is important is at the He/Fe boundary. 

We next solve for the shear rates by solving equation 
with either vq or v\, depending on which case is appropriate. 
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FIG. 2. — Accumulating models for accretion rates of O.lm^M (dotted 
line), 0.3niEdd (short-dashed line), and l.OflJEdd (solid line). For each model 
the base of the helium layer is taken to be at the unstable triple-« ignition 
depth where de^ a /dT = dt mo \/dT (thick long-dashed line in top panel). The 
top panel plots the temperature, and the solid circles mark the top and bottom 
of the convective zone. The middle panel plots the thermal diffusivity, K, 
and magnetic diffusivity, r). The bottom panel shows the Brunt-Vaisala fre- 
quency N for both thermal (Nt, lines) and compositional (JV„, solid squares) 
contributions. 



We assume angular momentum transport does not affect the 
thermal or compositional structure, in other words, that the 
transport is just happening "in the background." This allows 
us to assess which effects are crucial for subsequent iterations 
that include angular momentum transport in the actual struc- 
ture calculation. 

In Figure [3] lines denoted the profiles calculated assuming 
case 1 of the Tayler-Spruit dynamo. We assume that the con- 
vection is effectively instantaneous in transporting material 
and angular momentum, thus all the quantities become very 
small in the convective region. This is appropriate since the 
convective overturn timescale, H /Vconv ~ 10 s, is much less 
than the timescales for mixing or accretion. The squares in- 
dicate the corresponding values due to the compositional dis- 
continuity found by using case 0. Since the viscous timescale, 
fvisc = H 2 / v, is much less than the timescale it took to accrete 
to this base column f acc = y acQ /rh ~ 10 3 - 10 s s (for the range 
of m considered) our steady-state assumption is valid. 

These calculations highlight the importance of the composi- 
tional jump since f v j sc , q, and dE/dlny are all amplified here. 
This is because the large buoyancy reduces angular momen- 
tum transport across this boundary. The viscous time is nearly 
independent of accretion rate in regions where Nt 3> A^ M be- 
cause the viscosity in case 1 is independent of q. The energy 
deposition is always much smaller than the heat coming the 
crust, and it falls off somewhat faster with depth than expected 
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FIG. 3. — Angular momentum transport in the NS surface layers, for a spin 
of H = O.IQk, and accretion rates of O.lniEdd (dotted line), 0.3mEdd (dashed 
line), and l.OmEdd (solid line). The panels display (from top to bottom) the 
viscous timescale for angular momentum transport across a scale height, ? v ; sc , 
the shear rate, q, and the viscous energy deposition per logarithm column, 
dE/dlny. The squares indicate the corresponding values due to the compo- 
sitional discontinuity at the base of the accumulating layer. 



from the estimate presented in 33.31 due to electron degener- 
acy effects decreasing Nj. New heat sources at a depth of 
« 10 12 g cm" 2 could ease the difficulty calculations have in 
recreating the low igni tion columns needed to explain super- 
burst recurrence times dCumming et al.l20 06). but this viscous 
heating is not nearly enough to correct this problem. 

In Figure [4] we plot the magnetic fields found within the 
radiative zones for the models from Figure [3] We assume 
that within the convective zone the dynamo is not able to 
operate, and do not calculate a magnetic field here. It is 
also interesting to co mpare these fields to those derived by 
ICumming et al.l (1200 ll) . who calculated the steady-state fields 
expected when Ohmic diffusion balances advection through 
accretion. Their principal result was that the steady-state hori- 
zontal field drops by w m/0.02»?Edd orders of magnitude from 
the crust up through the ocean. In contrast, the magnetic 
fields we find are nearly constant with depth. We therefore 
do not expect these fields to persist if the accretion ceases for 
a time and instead to be expelled on an Ohmic diffusion time 
(~ days near the top of the ocean). When the dynamo is ac- 
tive, the steady-state is reached quickly enough that Ohmic 
diffusion can be ignored. 

In Figure|5]we compare what happens as the spin is changed 
by plotting n = 0.03, 0.1, and 0.3ft K (67, 220, and 670 Hz), 
all for m = O.lmEdd- Note that the shearing is most dramatic at 
smaller ft. In fact, the shearing and heating profiles are very 
sensitive to the value of f2, as was demonstrated by the ana- 
lytic estimates. We do not plot the associated magnetic field 
for these models since we have already plotted some exam- 
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FIG. 4. — Azimuthal and radial field components for a spin of f2 = O.IQk- 
Lines have the same meaning as in Fig. [3] 




Column Density (g cm 2 ) 

FIG. 5. — Same as Fig. [3] but for spins of 0.03Qk (solid line), 
(dotted line), and 0.3Ck (dashed line), all for m = O.lmEdd- 
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pies in Figure [4] and the analysis of 33.21 provides adequate 
estimates. 

These plots of q show that very little shearing is present. 
To emphasize this fact, in Figure [6] we plot the actual spin 
frequency as a function of depth found by integrating q for 
a range of accretion rates. The discontinuity at the He/Fe 
boundary complicates this estimate. Noting that 



n= / qSld\nr = - 



qtt—dlny, 



(54) 
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FIG. 6. — The spin frequency as a function of column depth for f2 = O.IQk, 
with m = 0.1 (dashed line), 0.3 (dotted line), and l.OntEdd (solid line). The 
dot-dashed line shows a constant spin frequency for comparison. 

we approximate the spin change at this boundary as 

AQ&qQH/R. (55) 

The change of spin is generally < 0.5 Hz across the accumu- 
lating layer, with the majority of the spin change occurring 
at the compositional boundary. The layer is very nearly in 
uniform rotation. 

4.2. Mixing and the Compositional Barrier 

The above calculations confirm that the viscosity is too 
large for either a pprec iable shearing or viscous heating. As 
demonstrated in 93.41 mixing should be important, but Fig- 
ures |3]and|5]argue that we must take into account the large 
at the He/Fe boundary. Assuming that scales like equation 
d53t . we estimate 

"'^(Aln^/O^) 1 / 2 . 



1.6 x 10 6 s _I fi\^r a 



(56) 



We substitute this into vq to solve for q using equation dU, 
which is used to estimate a mixing timescale at the boundary, 

t mix = H 2 /(a T sD ) 

/Aln^ I/3 



, n ,„3 -1 -5/3 4/3^5/3 .-4/3 

1.7 x 10 s a TS ^ L 33 p 6 T 8 m { fi . 



V 0.44 



(57) 



Using equation d48| > we find a ratio of 



1 -2/3 1/3^2/3 .-1/3 



; 0.77 a TS M L33 p 6 T 5 



'o.i 



Alnp 
0.44 



1/3 



(58) 

Unlike in equation d49l ), this new ratio depends on m. This is 
because the viscosity in case has a different dependence on 
q. Now as m is increased ? m ; x decreases faster than f acc . Above 
a critical accretion rate of 

'Aln^ x 



«i crit ,i =4.6 x 10 2 m Edd a T sML33P6 7^^0.1 



0.44 



(59) 

mixing can no longer be ignored. For densities and tempera- 
tures expected at the base of the accumulating layer, this crit- 
ical accretion rate lies in the range of 0. 1 - 1 .0»?Edd- Note that 
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FIG. 7. — Comparison of the mixing timescale, / m j x (solid line) versus the 
accretion timescale f acc (dashed line) for m — 0.1,0.3, and l.OmEdd (top to 
bottom panel; all with Q = O.ISIk and ots = D- The solid square denotes 
? m j x due to the compositional discontinuity. At low m, the squares are above 
the dashed line, demonstrating that ? m ; x > t RCC at the base of the accumulating 
layer, which prevents mixing to larger depths. At sufficiently high m, f m ; x < 
face at depths below the accumulation depth, so that mixing between helium 
and iron can occur. 

this depends very strongly on the spin rate, m cr i tj i oc f2 3 , thus 
we expect the slower spinning NSs to be considerably more 
affected by mixing. 

To test these analytic estimates we compare the mixing and 
accretion timescales in Figure [7] From the top panel to the 
bottom panel we increase m (fixing f2 = 0.1 J7k)- At low m, 
the He/Fe boundary at y acc (shown by the filled squares) acts 
as a barrier to mixing since f m ; x > f acc at this depth. When 
this happens it is a good approximation to ignore mixing and 
assume two separate layers during the accumulation phase. 
As m increases, f m ; x at y acc becomes less and less until fi- 
nally f m ; x < f acc , so that material should be mixed past y acc . 
When this occurs, our accumulating model can no longer ig- 
nore mixing. The mixing between helium and iron occurs 
down to a depth where f m ; x is equal to the length of time ac- 
cretion has been taking place, f acc = y^jm. The key point we 
want to emphasize is that because of the buoyancy barrier, the 
effect of mixing turns on abruptly, and when it does, mixing 
will occur well past y. dcc . 

5. THE EFFECTS OF TURBULENT MIXING 

Once f m ; x < f acc at the buoyancy barrier, the compositional 
profile of the NS is very different, which we now explore. 
We treat the mixing as complete, which we diagram in Fig- 
ure [8] and summarize here. Material accretes at a rate m for 
a time f acc with a helium mass fraction Yq, supplying a col- 
umn of material y acc = mt- dcc . The total column of helium that 
accretes during this time is therefore Yoy acc . Mixing causes 
this newly accreted material to mix past y acc to the mixing 
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FIG. 8. — Diagram demonstrating the main features of turbulent mixing. 
Material mixes down to where f m j x = f acc , which defines the depth y m j x . The 
total amount of material that accretes is y aC c = «i/ aC c , giving a total accreted 
column of helium of Foyacc . Within the mixed layer the helium mass fraction 
is diluted to new mass fraction of Y mlK = Jo^acc/ymix- 



depth y m [ x defined as where f m ; x = f acc . The helium is fully 
mixed down to this depth, resulting in a diluted mass fraction 
JWx = Yoy acc /y m i x within the mixed layer. In the following 
sections we consider two scenarios that can result from the 
mixing: (1) unstable ignition of the mixed fuel, and (2) stable 
burning when the material mixes to sufficient depths. 

5.1. Numerical Calculations of Mixed Ignition 

As the column of accreted material grows, it can still reach 
the correct conditions for unstable ignition, but mixing causes 
two changes: (1) the ignited layer has a diluted helium frac- 
tion, Ymix < Yq, an d (2) the ignition takes place at the base 
of the mixed layer, at a depth y m ; x > y acc , resulting in a re- 
currence time for ignition much less than when mixing is not 
included. 

Both of these effects are easiest to explore using a semi- 
analytic model. In this section all calculations use Yq = 1 . We 
consider other values of To for our analytic estimates in the 
next section. We solve for the mixed accumulating structure 
by first assuming an amount of accretion y acc , which for a 
given m implies an accretion time y- da: /m. We integrate the 
radiative diffusion equation (eq. J52]) down to a depth where 
fmix = face giving y m i x , where we assume a constant flux pro- 
file as was discussed in $4] since little helium burning takes 
place during accumulation. We then estimate the mixed he- 
lium fraction as Y m i x =y acc /y m i x . This estimate is improved by 
substituting Y m [ x back into our envelope integration, and iter- 
ating until T m i x converges. The non-helium component in the 
layer is taken to be iron. 

In Figure [9] we plot the resulting profiles for a NS accreting 
at m = O.lmEdd and ft = O.IS^k- In the four panels we consider 
values of y acc of 10 6 , 3 x 10 6 , 10 7 , and 4 x 10 7 g cm" 2 (from 
left to right and then up to down, denoted by the filled circles), 
which is meant to mimic the accumulation of fuel on the NS 
surface. For each integration the envelope profile continues 
down to a depth ymi x , which also gives F m i x as displayed in 
each panel. The last model reaches the conditions necessary 
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FIG. 9. — The four panels show how the fully mixed accumulating layer 
evolves in time until it reaches conditions necessary for unstable ignition. 
The parameters of the NS are m = O.lmridd an d = O.IQk- In each panel, 
the column of helium that has been accreted is denoted by a filled circle, 
which is from left to right, and up to down y acc = 10 6 , 3 X 10 6 , 10 7 , and 
4 X 10 7 g cirr 2 . Mixing takes place down to the column reached by the thin 
solid line. The mixed helium fraction, F m j x is displayed in the upper left-hand 
corner of each panel. The ignition curve associated with each F ra j x is shown 
as a thick dashed line. 



for ignition at the base of the mixed layer. The recurrence time 
for these mixed-ignition models is much less than for those 
without mixing. For the plotted model, the recurrence time is 
f rec = 4 x 10 7 g cm~ 2 /1.5 x 10 3 g cm" 2 s" 1 « 44 minutes. In 
contrast the model shown in Figure|2]with the same accretion 
rate of m = Q.lniEdd has f rec w 1.5 days! The shorter recur- 
rence time is not only due to mixing carrying helium down to 
deeper depth, but also a change in the thermal profile. Since 
significant iron is mixed up into the accumulating material, 
the free-free opacity, which scales Kff cx Z 2 /A, where Z and 
A are the charge and mass per nucleon, respectively, is now 
the dominate opacity mechanism. The accumulating layer is 
more opaque and therefore hotter for a given flux in compari- 
son to the pure-helium models considered before, which con- 
tributes to the shallow ignition depths. 

We calculate f rec for a grid of models with various m and 
SI in Figure [10] The recurrence time is shorter for stronger 
mixing, which occurs at high m or low SI, and can be as short 
as « 5 - 30 minutes. We warn though that all of these cal- 
culations assume that complete mixing can occur, and as we 
already showed in $4] buoyancy may prevent this (eq. f59l ). 
Nevertheless, it is interesting to calculate the mixed-ignition 
conditions for a wide range of parameter space because the 
conditions left from previous bursts may vary, and quantities 
such as A In /.i may be smaller at times if, for example, there is 
incomplete burning in a previous burst. At sufficiently high m 
or low SI the ignition takes place at high enough temperatures 
that the envelope does not ignite unstably. We consider this 
case in more depth in the following section. 

The short recurrence tim es that we find are similar to some 
seen for multiple bursts (Gall oway et alj|2006l and references 
therein). However, our model cannot explain the energetics 
of these bursts. The energetics are typically quantified in a 
distance independent measure, the so-called a-value, which 
is the ratio of energy released in the persistent emission be- 
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FIG. 10. — The recun'ence time for mixed-ignition models as a function of 
m. The symbols denoted different spins, as shown in the key. Models that 
are at sufficiently high m or low Q do not ignite unstably, and thus are not 
plotted. 

tween bursts to the energy of the burst itself. For pure helium 
ignition a ~ 100. In contrast, the q- values for the shor t recur- 
rence time bursts are typically ~ 10 (Boi rin et al.l20 07). Since 
the only helium that will burn in our mixed ignition models is 
that accreted since the last outburst, we always find a ~ 100. 
Therefore we still require an additional nuclear energy source, 
such as incomplete burning from the previous X-ray burst, to 
explain such a low a-value. This in fact may not be a problem 
because incomplete burning may naturally explain the short 
recurrence times, since this would lead to smaller composi- 
tional gradients and therefore stronger mixing. 

We plot the Y m [ x as a function of m and f2 in FigureQl] The 
crosses correspond to models that ignite unstably. The filled 
circles are stable accreting models that are discussed on 35 .31 
and 35.41 The implications of this mixed ignition for the un- 
stable burning during the X-ray burst can be easily tested by 
more sophisticated numerical simulations by just considering 
a mixed accumulating column, for example, by artificially set- 
ting Y w 0.1 -0.6 as shown in this Figure [TT] Since viscous 
heating is negligible, these initial tests do not need to resolve 
the shearing profiles. 

5.2. Analytic Estimates of Mixed Ignition 

We now estimate the properties of the mixed ignition mod- 
els. These solutions directly show how mixing depends on the 
properties of the accumulating layer, in particular the pref- 
actor cuts and Yq, without having to consider a multitude of 
models. The effects of a free-free opacity are important in 
deriving the correct atmospheric conditions. To include this 
analyt ically, we use the free-free opacity from Scha tz et"aT] 
(1999), simplified to a one-component plasma and with the 
dimensionless Gaunt factor set to unity, 



_ __ 2-1 A> 
Kff«3.77cm g 



T V 2 A 

i 8 



(60) 



where we have estimated p e ps 2 as is correct within 10% for 
the any of the elements of interest. Integrating the radiative 
diffusion equation (eq. 1521 ). assuming a constant flux of F = 



10 21 ergs cm 



s F21 and an ideal gas equation of state, the 



FIG. 11. — The mixed helium mass fraction, y m i x as a function of m and Q. 
The crosses are mixed ignition models. The filled circles are for stably burn- 
ing models and give the amount of helium present in the steady-state mixing 
and burning layer. Steady-state burning therefore require either high m or 



low Q. The scaling of Y m 



is derived in equation (69) and is consis- 



tent with these numerical results. Without mixing all of the these considered 
models ignite unstably since stable accretion requires m > lOniEdd- 



temperature as a function of column y is 

7\y)=1.8xl0 8 K (mi.33F 21 Z 2 /A) 2/17 - 4/17 



- V 8 



(61) 



where = y/10 8 g cm 2 . Using equation ((47), we find the 
mixing timescale as a function of y, 



: 3.5 x 10 3 s c&tfg (FnZyAr -m^yV 



(62) 



Note that f m j x cx y is a higher power than f acc =y/m (xy. 
This explicitly shows that mixing dominates at lower y, but 
accretion always wins at some depth. Setting f m ; x equal to 
f acc = Jacc/'w gives the depth where mixing can extend to for a 
given column of accreted material y acc , 



Jmix0'acc)= 1-6 X 10 8 g CHI ^tsVlM 



x (^zva) - 1 Vi 55 y 



0.55„0.73 

acc,8' 



(63) 



where y acc .8 = jacc/10 8 g cm 1 . The mixed helium fraction 
down to this depth is 



^mix^acc) 



lay a, 



3 7 mix(3 7 acc) 

-0.63 a TS roMi.33 (^21^ I A) S2 .i ^acc.s- 



(64) 

Since Y m i x cx y ac 27 the strength of mixing decreases (Y m \ x gets 
larger) as more material accreted, which was demonstrated by 
the four panels in Figure [9] 

We next estimate what ignition depth is expected for this 
fully mixed accumulating layer. The energy generation rate 
for triple-a burning is approximated as 



n 2 Y 3 

e 3a = 5.3 x 10 23 ergs g" 1 s" 1 



exp 



44 



, (65) 



where / is factor that accounts for screening effects. To 
make progress analytically we expand the exponential as 
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eX p(-44/r 8 ) w 7.95 x 10" 10 (r 8 /2.1) 21 . Using our tempera- 
ture profile (eq. lloTl ), the condition that de3 a /dT = de coo \/dT 
(eq. l50l ) implies an ignition depth of 



3V> : 



= 9.4xl0 7 gcm- 2 r 0I VL33W 13 



' 21 



x(z 2 /A) -0.28 F -0.44^ 



(66) 

Setting y lgn = y m i x from equation ( |63t , we solve for y acc , the 
critical column of material that must be accreted to cause ig- 
nition. This is then substituted back into equations d63l and 
d64l i to find that ignition occurs at a depth 
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(68) 



Comparing this with the numerical calculation is easiest if we 
assume F is set by m (F21 = 2.2;«o.i) as well as scaling Z 2 /A w 
12 and /x w 2.1 as appropriate for the iron-rich composition. 
This gives 



, n on -0-86 --0.21n0.65 



(69) 



The recurrence time is f lec = T m j x .i gn y m j x .ig n /m, resulting in 
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Equations d69l and (TTOb confirm the scalings found for the 
numerical calculations in Figures [10] and [TT] We have plotted 
the F m i x oc m~° 21 scaling in FigureQT]to emphasize this. These 
analytic results also show how strongly these results depend 
on the parameter axs, for which Tmi^ign is especially sensitive. 

5.3. Steady-State Mixing and Burning 

Another possibility is that the helium is mixed and burned 
by triple-a reactions in steady-state, leading to stable burning. 
The basic idea is similar to that described above for mixed 
ignition, except now the depth of the accumulating layer is 
set by the helium burning timescale, ha = ET, a /(YeT, a ), where 
Ej a = 5.84 x 10 17 ergs g" 1 is the energy per mass released 
from this burning. As shown in Figure [12] material is mixed 
to sufficient depths where ha (which decreases with depth) 
is equal to f m j x (which increases with depth), which defines 
a mixing (or burning) depth y mix = y(f mix = t ia ). During a 
mixing timescale, the amount of material that is accreted is 
y acc = mfmix^mix), so that the total column of helium that has 
been accreted is 5 / oJ7if m i x (y lmx ). This is diluted over a depth 
y m ; x , so that the mixed helium fraction is y m i x = ioJacc/jmix- 
This is all occurring in steady-state, material moves through 
the mixed layer at a rate m, but this layer does not move up 
or down in pressure (column) coordinates, as the burning is 
stable. 

The basic results of steady-state mixing and burning are 
best shown using a simple numerical model. The equation 
that describes helium continuity, including depletion by triple- 
a bu rning, becomes, in the plane-parallel limit dFuiimotol 
1993), 
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We assume changes in p and D with depth are small in com- 
parison to changes in Y, and, following our derivation of the 
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FIG. 12. — Diagram demonstrating the main features of steady-state mix- 
ing and burning. Material can mix further down to where t mix = t^ a , which 
defines the depth y m i x . During a timescale f ra ix(}'mix) the amount of material 
which has accreted is Vacc = "jfmix(ymix)- Within the mixed layer the helium 
mass fraction is diluted to new mass fraction of Y m j x = Yoyace /ymix- 



angular momentum equation in 
limit to derive 



mdY _ D d 2 Y 

p dz dz 2 



we take the steady state 



Y 

ha 



(72) 



In the limit where f acc ^> f m ; x that we are interested in, the 
term on the left-hand side can be dropped. Finally, making 
the approximation that pDd/dz ~ y/t mix, we find 



dY 
dy~ 



ha y 



(73) 



This equation mimics the properties we expect from mixing. 
When mixing is strong f m ; x / h a <C 1 , and dY jdy rj 0, the com- 
position does not change with depth. At the depth where 
fmix ~ hat dY jdy ~ —Y/y and the helium is depleted expo- 
nentially. All the helium burns into carbon, so that carbon has 
a mass fraction X\% = l-Y. 

The envelope profiles are found from simultaneously in- 
tegrating three differential equations: (1) radiative transfer, 
equation d52"l l. (2) the entropy equation, dF/dy = -£3 Q , and 
(3) continuity of helium, equation ( l73l . We integrate using 
a shooting method, but first we must set three boundary con- 
ditions for the flux, temperature, and helium mass fraction. 
Since all of the accreted helium must burn if the envelope is 
in steady-state, we set the surface flux to F = Yp Ei, a m + F c , 
where F c = 150 keVnuc _1 (m) (as discussed in 134.11 ). The 
surface temperature i s set from the radiative zero solution 
dSchwa rzschild 1958). The helium abundance in the mixed 
region, F m j x , is an eigenvalue. It is varied until shooting gives 
the correct base flux of F c . This is easily found through itera- 
tion since when F m ; x is set too large, too much burning occurs, 
and the base flux is too small (and vice versa for small Fmi x ). 

In Figure [13] we plot a steady-state envelope using Yq = 1, 
m = 0.3f7iEdd, fi = O-lf^K, and axs = 1. The shooting method 
demonstrates that the initial helium abundance within the 
mixed layer is Y m [ x = 0.097. The top panel shows the tem- 
perature profile (thin solid line), and the critical curve for sta- 
bility where de^a/dT = de coo \/dT (thick dashed line). The 
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FIG. 13. — An example steady-state mixing and burning envelope model 
using Fq = 1, m = l.OnirJdd. = 0.03f2K, and a-rs = !• The material within the 
mixed layer has Y ma = 0.097. The top panel shows the temperature profile 
(solid line), as well as the critical curve ignition using a helium mass fraction 
of 0.052 (thick dashed line, which is the helium mass fraction at the burning 
depth). The middle panel compares the key timescales. The bottom panel 
shows the helium (dashed line) and carbon abundances (dotted line), as well 
as the energy generation rate for helium burning (thick solid line). 



bottom panel shows that the majority of triple-a burning oc- 
curs at « 3 x 10 s g cm" 2 (thick solid line), which is where the 
helium is depleted. Comparing the middle and bottom panel 
shows that the majority of the burning takes place near where 
ttsux = ha (as required by construction), which is deeper than 
where f acc = tj a (the normal condition for steady burning). 

To make sure the steady-state models we find are physically 
realizable, we must check the thermal stability of the helium 
burning. In Figure [l4]we compare the quantities de-3 a /dkiT 
and afe coo i/iilnr for three different accretion rates. If the 
cooling derivative is always larger, then the model is thermally 
stable. We find stable accretion at m's considerably less than 
the stable accretion rate of m « lOmEdd ~ 2 x 10 6 g cm -2 s -1 
expected for pure hel ium accretion estimated without mixing 
jBildstenll 1 995L 1 1 998a|) . Models that are found to be stable in 
this way are plotted as filled circles in Figure QT| (from 35. 11 1. 
Only at large m and small 57 are the models found to be stable. 
If we where to increase aj$, a wider range of the models in 
Figure QT]become stable. 

5.4. The Critical rhfor Stability 

Since we have found a set of models that can stably accrete, 
mix, and then burn helium, it is interesting to ask what accre- 
tion rates and spins are required for this to occur, and how 
does it depend on parameters such as ar$- 

First we must derive the correct condition for stability in- 
cluding the fact that free-free opacity is important in setting 
the radiative profile. The one-zone condition for stability at 
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FIG. 14. — A comparison of dej a /dlnT (solid lines) and de coo i/dlnT 
(dashed lines) as a function of depth for m = 0.1, 0.3, and l.Omrjdd (top to 
bottom panel). All the models use f2 = 0.03f!K, «ts = 1> an d 5o = !• This 
demonstrates that only the m = l.OniEM model is stable out of these three. 



the base of the mixed layer is 



d€3 a da c ool 

dlnT < dWT 1 



(74) 



where the derivatives are taken at constant pre ssure. For e^g o c 
T^p x , where ( = 44/^-3, stability requires (Bild stenll998al) 
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(75) 



Substituting the scalings for a free-free opacity and ideal gas 
equation of state we find 7g > 3.26 is required for stability. 

By substituting the analytic form we found for the mixed 
ignition depth (eq. l67l ) into the temperature profile (eq. l6Ti ) 
we can find the temperature at the base of the mixed layer, 



r ign = 2.54xl0 8 K4 s °Vi!^o i 067 



(76) 



By simply asking when 7i g 
bilizing m of 



> 3.26 x 10 8 K, we derive a sta- 



mcrit,2=1.0m Edd a^ s 83 ^;f . (77) 

This is in reasonable agreement with Figure [TT] which shows 
stability occurs for 0.9mEdd <th< l.(Medd for f2 = 0.1 Ok- It 
is interesting that »i cr ; t 2 is near (within an order of magnitude) 
a value where bursts are observed to change. It is conceivable 
that this mechanism may act to stabilize X-ray bursts for the 
hydrogen-rich accreting systems, and given the strong scaling 
fWcrit,2 has with cuts it is possible that more detailed calcula- 
tions could give results that agree even better with the critical 
accretion rates that are observed. We discuss this idea in more 
detail in the following section. 
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6. DISCUSSION AND CONCLUSION 

We have revisited the problem of angular momentum trans- 
port in the surface layers of accreti ng NSs. We foun d that the 
hydrodynamic instabilities used by Fuii motol d 19931) in a pre- 
vious study are dwarfed by the magnetic effects of the Tayler- 
Spruit dynamo. The large viscosity provided by this process 
results in a very small shear rate and negligible viscous heat- 
ing. The turbulent mixing is sufficiently large to have impor- 
tant consequences for X-ray bursts. We constructed simple 
models, both analytic and numerical, to explore mixing for 
pure helium accretion. From these models we can make a few 
conclusions that are likely general enough to apply to most 
viscous mechanisms. As a guide, we show the different burn- 
ing regimes we find in Figure [15] These can be summarized 
as follows: 

• Mixing is strongest at large m (when angular momen- 
tum is being added at greater rates) and small SI (which 
gives a larger relative angular momentum between the 
NS and accreted material). 

• Mixing has trouble overcoming the buoyancy barrier 
at chemical discontinuities. But once mixing breaks 
through this, it extends down to a depth where f m i x = 4cc, 
which is generally much deeper than the accreted col- 
umn. This means that mixing should turn on abruptly 
(as a function of either m or Jl). It also means that the 
importance of mixing depends on the particular ashes 
left over from previous bursts. If, for example, incom- 
plete burning results in small compositional gradients, 
mixing would be important in subsequent bursts. 

• Mixing of freshly accreted material with the ashes from 
previous X-ray bursts can lead to two new effects. First, 
the layer may ignite, but now in a mixed environment 
with a short recurrence time of ^ 5 - 30 minutes. Sec- 
ond, if the mixing is strong enough, accreted helium 
can mix and burn in steady-state, quenching X-ray 
bursts. Both of these regimes have observed analogs, 
na mely the short rec urrence time bursts (for example 
in, iBoirin et~aT]|2007l) and the stabiliza tion of bursting 
seen at w 0.1/«Edd dCornelisse et alj|2003l) . 

The mixed ignition case can be studied easily using the cur- 
rent numerical experiments (e.g JWooslev et al. 2004) by just 
artificially accreting fuel with a mixed composition of w 
0.1-0.6. These calculations are simplified by our conclusion 
that shearing and heating can be ignored, at least for initial 
studies. We next conclude by speculating about some of the 
other ramifications of turbulent mixing. 

6.1. Superbursts 

An ongoing mystery in the study of bursting NSs is 
the recurrence times for superbursts, thermonuclear ig- 
nition of carbon in the X- ray burst ashes at c o lumns 
of sa 10 11 - 10 12 gem" 2 dCumming & Bildstenl 120011: 
IStrohmaver & Brownl 12002b . This problem could be alle- 
viated by enha nced heating from th e core on the order of 
1 MeV nuc" 1 (Cumming et al. 2006), but this is more than 
is expected theo retic ally, even in the newest calculations 
(Gupt a~et al.l2006l) . Shear heating is not large enough to solve 
this problem, as demonstrated in Figures [3] and [5] 

The regime of stable helium burning we have found may, 
however, create a carbon rich ocean that would assist in the 
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FIG. 15. — Summary of the three regimes of burning found for models 
including mixing. The boundaries between each regime are shown as shaded 
regions to emphasize possible uncertainty in the strength of mixing (we con- 
sider «xs = 0.7 — 1.5). The heavy shad ed region divides where the buoyancy 
barrier is overcome (m cr ; t i, eq. 1591 ), and the light sh ade d region divides 
between stable and unstable mixed burning (m cr jt,2, eq. 1771 ). An additional 
uncertainty in m alt \ is the value of A In fi, which could vary depending on 
the results of previous bursts. 

ignition of superbursts. Carbon fractions of greater than 
10% are needed to re produce the lightc urves and recurrence 
times of superbursts ( Ke ek et all 12006). Calculations of rp- 
process burning show th at unstable burnin g cannot give car- 
bon fractions this hig h (ISchatz et alJl2003l) . Observat iona lly, 
lin'tZand et al.l d2003l) showed that the a-value (see jj5.U is 
preferentially large for superbursting systems, indicating that 
some stable burning is occurring, perhaps due to the turbulent 
mixing we have studied. 

6.2. Hydrogen-rich Accretion 

One of the main deficiencies of our calculations are the sim- 
plified compositions, since most accreting NSs are expected 
to be accreting a fuel abundant in hydrogen. If the only cru- 
cial burning is triple-a, such envelopes can be considered 
within the framework of our models by using a solar value 
of Yq = 0.3. Our models fail, though, when sufficient carbon 
is produced by triple-a to feedback into hydrogen burnin g 
(which burns via the hot-CNO cycle. iHoyle & Fowlerll 19651) . 
Such a scenario is interesting thought because it could poten- 
tially produce very hydrogen poor bursts. 

In the standard theoretical framework for X-ray bursts, 
flashes should be hydrogen poor at low m (when there is suf- 
ficient time for the hot-CNO cycle to act), and then become 
mixed hydrogen-helium fuel at higher m. Paradoxically, ob- 
servations show a transition at a seemingly universal luminos- 
ity of » 2 x 10 3 7 ergs s" 1 (approximate ly O.lmEddX but in the 
opposit e sense dCornelisse et alj [20031) . ICooper & Naravanl 
(2006b) argue that their models in fact show this transi- 
tion (also see [Narayan & Heyll 120031) because carbon cre- 
ated during helium simmering increases the the hot-CNO 
rate and the temperature, which decreases the temperature 
sensitivity of triple-a reactions. To explain the discrep- 
ancies betw een their models and more detailed n umerical 
simulations (Woosley et al. 2004; Hegeret al. 2005) requires 
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a decrease in the breakout reactions rate of 1 5 Q(a,7) 19 Ne 
(ICooper & Narayanl l2006at iFisker et al l 120061) . Unfortu- 
nately, the most recent experimental results do no t support 
a decreased rate dTan et al1l2007b IFisker et al. 2007). Further- 
more, helium accreting systems show this same transition in 
bursting properties, which is not explained within their frame- 
work. 

Another idea that may recreate this trend is that the frac- 
tion of the star covered by the fuel increases with the global 
accretio n fast enough th at the local accretion rate actually de- 
creases (Bildsten 20 00|). This interpre tation is supported by 
the recent work of iHeger et ail (1200 5!). who claim that the 
mHz oscillations observed at a round this same univer sal lu- 
minosity of 2 x 10 37 ergs s" 1 (iRevnivtsev et al.ll2001l) are a 
probe of the local accretion rate where bursting is occurring. 
The main problem with this suggestion is that it is difficult 
to understand how the covering of the surface can remain so 
anisotro pic all the way down t o the d epths of where ignition 
occurs. Inogamov & Sunyaevl (1 19991) calculated the spread- 
ing of accreted material from the equator, where the viscosity 
is due to a turbulent boundary layer, and find that spreading 
occurs orders of magnitude shallower in depth than where ig- 
nition takes place. 

Since mixing gets stronger with m, it may be that the ob- 
served transition is in fact hydrogen being turbulently mixed 
and burned analogous to what we have found for pure- 
helium accretion. A possible implication of such an expla- 
nation is that slowly spinning NSs are more likely to have 
their hydrogen depleted resulting in helium-rich bursts. As- 
suming that the burst oscillation frequencies are indicative 
of the NS spin frequency (which i s close to true or ex- 
actly true for all current explanations lHevlll2004t lLedl2004t 
iPiro & Bildstenl 120051: iPayne & Melatosl l2006l) . the systems 



can be broken into slow spinning (~ 300 Hz) and fast spin- 
ning (~ 600 Hz) classes. This slowly spinning clas s in- 
cludes 4U 1916- 053 (270 Hz. iGaliowav et all bOQlt). 4U 
1702-4 29 (330 Hz.lMarkwardt et all 1999b and 4U 1728-34 
(363 Hz. IStrohmayer & Markwardtll 19991) . Do these systems 
show helium-rich looking bursts, and if so, is this due to turbu- 
lent mixing destroying the hydrogen they are accreting? This 
can be answered by l ooking at the re c ent su mmary of RXTE 
burst observations by iGalloway et all (120061) . 4U 1916-053 
has bursts consistent with helium-rich fuel, but this system 
also has an orbital period of rj 50 min dGrindlav et al.l ll988). 
This is an "ultracompact" system in which the donor is too 
small to be a H-rich star, so the accretion is probably helium- 
rich to begin with. The other two NSs both have bursts with 
decay times and a-values that suggest helium-rich fuel, which 
is also supported by model fits to radius expansion bursts 
from 4U 1728-34 (faallowav etaill2006l) . Furthermore, the 
bursts of 4 U 1728 — 34 have look very similar to those of 4U 
1820-30 je umming 2003), a know n ultracompact (see dis- 
cussion in Podsiadlowski et al. 2002, and references therein). 
While it is possible that 4U 1702-429 and 4U 1728-34 
have helium-rich donors (and are thus ultracompacts), it may 
also be that they are accreting hydrogen-rich fuel and show 
helium-rich bursts because their low spins lead to mixing. As 
binary parameters of these systems become better known, it 
will be more clear whether turbulent mixing is indeed needed 
to explain their burst properties. 
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APPENDIX 

MIXING AND ENERGY CONSERVATION FOR THE TAYLER- SPRUIT DYNAMO 

In jQ.ll we summarized the prescriptions that lSpruitl d2002l) provides for the Tayler-Spruit dynamo. The mixing diffusivity is 
assumed to be equal to the turbulent magnetic diffusivity. Although this is plausible, it is not shown rigorously. Below we argue 
that such a scaling is consistent with energy conservation. 

Consider a layer differentially rotating with a speed AV. The equation for the energy per unit mass expresses the fact that 
energy from differential rotation can go into either viscous shearing or vertical mixing 



d 
dt 



(Avy 



1 



2ApgH 

2 oi\r2; 



where D = H 2 /t n 
1988) 



-2N Z D, (Al) 
is the mixing diffusivity. The ratio of the two right hand terms is called the flux Richardson number (Fujimoto 



4N 2 D AD 

R f=^T^2 = — RL (A2) 
i/(qliy v 

where Ri = N 2 /(q 2 £l 2 ) is the Richardson number. The flux Richardson number is interpreted as the ratio of energy that goes into 
mixing versus heating (and usually taken to be Rf ~ 0. 1 — 1). From this result we can solve for the mixing diffusivity 



D = v. 

4Ri 



(A3) 



This matches (up to a factor of order unity, Rf/4-) what Spruit (2002) gives for the mixing diffusivity. The factor of Ri in the 
denominator takes into account the difficulty in overcoming buoyancy to do mixing. The larger the buoyancy, the larger Ri is and 
the smaller the mixing diffusivity. 
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